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ABSTRACT 


This  Semiannual  Technical  Summary  covers  the  period  1  October  1980 
through  31  March  1981.  It  describes  the  significant  results  of  the  Lin¬ 
coln  Laboratory  Multi -Dimensional  Signal-Processing  Research  Pro¬ 
gram,  sponsored  by  the  Rome  Air  Development  Center,  in  the  areas  of 
image  segmentation  and  classification,  adaptive  contrast  enhancement, 
and  iterative  implementations  for  multi -dimensional  digital  filters. 
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MULTI -DIMENSIONAL  SIG  N  A L -  PR O C E  S SI NG 
RESEARCH  PROGRAM 


i.  INTRODUCTION  AND  SUMMARY 

The  Lincoln  Laboratory  Multi -Dimensional  Signal-Processing  Research 
Program  was  initiated  in  FY  SO  as  a  research  effort  directed  toward  the  de¬ 
velopment  and  understanding  of  the  theory  of  digital  processing  of  multi¬ 
dimensional  signals  and  its  application  to  real-time  image  processing  and 
analysis.  A  specific  long-range  application  is  the  automated  processing  of 
aerial  reconnaissance  imagery.  Current  research  projects  which  support  this 
long-range  goal  are  image  modeling  for  segmentation  and  classification,  tech¬ 
niques  for  adaptive  contrast  enhancement,  iterative  implementation  of  multi¬ 
dimensional  digital  filters,  and  multiprocessor  architectures  for  Implementing 
image  processing  algorithms.  Results  in  these  research  areas  over  the  past 
six  months  are  described  in  this  Semiannual  Technical  Summary. 

~s>  in  the  area  of  image  segmentation  and  classification,  we  have  been  devel¬ 
oping  a  hierarchical  segmentation  scheme  for  processing  images  with  several 
region  classes.  This  approach  appears  to  offer  improvements  over  a  direct 
multiclasa  segmentation.  Examples  are  given  in  Sec.  Z, 

^Adaptive  contrast  enhancement  techniques  have  proved  useful  in  several 
areas.  We  have  used  adaptive  contrast  enhancement  as  a  preprocessor  for 
image  segmentation  based  on  texture  rather  than  gray  level.  We  have  also 
applied  these  techniques  to  aerial  images  degraded  by  light  cloud  cover  and 
haze.  The  primary  effects  of  this  type  of  degradation  are  a  reduction  In  con¬ 
trast  and  an  increase  in  the  local  average  intensity  of  the  image. »  In  Sec.  3, 
we  show  examples  of  the  improvement  which  the  adaptive  contrast  enhancement 
techniques  afford  on  images  suffering  from  this  type  of  degradation. 

In  Secs.  4  and  5>  we]  discuss  tw*<  aspects  of  the  iterative  Implementation 
of  multi-dimensional  digital  filters,  mia  Hrat  concerns  cpatial  truncation 
effects  which  occur  during  the  iteration  because  the  image  frame  buffer  has  a 
finite  storage  capability.  The  errors  caused  by  this  truncation  are  closely  - 
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related  to  the  solution  of  a  boundary  value  problem  with  boundary  conditions 
specified  on  the  frame  edges.  These  errors  can  be  eliminated  by  including 
the  boundary  conditions  in  the  spatially  truncated  iteration. 

Section  5  discussed  potential  multiprocessor  architectures,  for  realizing 
multi-dimensional  signal-processing  operations  such  as  those  needed  in  the 
iterative  implementation.  An  underlying  problem  is  the  efficient  partitioning 
of  multi -dimensional  signal-processing  problems  among  the  several  proces¬ 
sors  in  a  multiprocessor  architecture.  e — 

2.  IMAGE  SEGMENTATION 

In  the  previous  Semiannual  Technical  Summary  Report,1  we  described  a 
segmentation  procedure  based  on  linear  filtering  methods  to  model  texture  in 
local  regions  of  an  image  and  a  Markov  random  field  to  model  region  transi¬ 
tions  within  the  image.  We  further  showed  application  of  this  method  to  Images 
containing  two  region  types.  The  segmentation  algorithm  described  in  Ref.  i 
ia  also  applicable  to  cases  where  there  are  more  than  two  types  of  regions  in 
the  image.  The  conditions  for  the  segmentation  are  as  follows. 

Assign  a  pixel  with  coordinates  (n,  ra)  to  class  1  where  i  is  the  class  for 
which 

E.  (n,ro)  ? 

-  +  in^  -21nPr[i|Sn  )  (i) 

ai 

is  minimum.  In  the  above  expression,  E,  is  the  error  in  linear  prediction  of 

*  2 

a  pixel  from  a  surrounding  set  of  pixels,  <r,  is  the  variance  of  the  prediction 
error,  and  Pr  { i  j  Sfj  ^  j  is  the  probability  that  the  given  pixel  belongs  to  class  i 
given  the  classes  of  a  surrounding  set  of  pixels.  (See  Refs.  1  and  2  for  more 
details.)  During  the  current  reporting  period,  we  have  been  applying  the  multi- 
class  form  of  the  algorithm  to  aerial  photographic  data  and  experimenting  v.ith 
strategies  for  the  segmentation  of  multiclasa  images.  Our  results  indicate 
that  it  is  often  better  not  to  attempt  to  segment  an  image  into  a  large  number 
of  classes  all  at  once.  Instead,  one  can  first  segment  the  image  into  a  few 
categories  and  then  perform  additional  segmentation  within  each  category.  We 


refer  to  this  approach  as  a  layered  segmentation  strategy.  Several  examples 
of  using  this  layered  approach  are  described  in  the  following  paragraphs. 

In  Ref.  1  we  stated  how  an  image  which  had  been  segmented  into  tree  and 
field  regions  could  be  further  segmented  into  regions  containing  only  large  or 
only  small  trees.  Figure  1  compares  results  of  this  layered  segmentation 
approach  to  a  nonlayered  approach  for  this  image.  The  image  in  Fig.  1  (left) 
was  first  segmented  into  the  broad  categories  of  trees  and  fields.  Within  the 
tree  category  we  further  segmented  the  image  into  large  and  small  tree  re¬ 
gions;  in  the  field  category  we  segmented  the  image  into  two  different  field 
types.  The  region  boundaries  arising  from  this  layered  approach  are  shown 
in  color.  (Note  that  the  small  field  in  the  lower  left-hand  corner  has  the  same 
texture  as  the  field  adjacent  to  it  and  thus  is  not  recognized  as  a  separate  re¬ 
gion  by  the  segmentation  algorithm.)  Figure  1  (right)  shows  the  results  of 
segmenting  the  image  into  the  four  regions  in  a  single  step.  The  regions  iden¬ 
tified  by  the  algorithm  are  coded  in  different  colors.  Figure  i  shows  that  the 
layered  segmentation  strategy  produced  a  more  accurate  result.  Although  both 
approaches  gave  quite  accurate  estimates  of  the  boundary  between  the  two  field 
types  and  the  boundary  between  the  trees  and  field,  the  layered  approach  was 
able  to  make  more  subtle  distinctions  between  the  small  and  large  trees.  Both 
segmentation  strategies  had  some  problems  near  the  boundaries  of  the  image 
due  to  the  discontinuities  arising  there  but  the  problems  were  more  apparent 
in  the  spot  near  the  left  edge  of  the  image  in  Fig.  i  (right).  The  results  of  this 
comparison  are  not  realty  surprising  since  in  a  four-category  segmentation 
there  is  more  room  for  error  in  the  classification  of  the  individual  pixels  than 
in  successive  two-category  segmentations.  A  few  points  in  the  small  tree  re¬ 
gion  classified  as  "field"  during  the  initial  classification  phase  of  the  algorithm 
can  propagate  their  effects  through  the  iteration  employed  in  the  algorithm  and 
lead  to  a  final  result  that  is  incorrect.  On  the  other  hand,  a  layered  approach 
may  produce  a  large  error  in  the  initial  categories  that  cannot  be  corrected  at 
successive  levels  of  segmentation.  Thus,  one  can  expect  that  a  layered  ap¬ 
proach  will  produce  best  results  where  there  is  clear  separation  between  the 
initial  categories  and  subtle  distinctions  between  the  classes  within  those  initial 
categories. 
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Fig.  i.  Comparison  of  a  layered  segmentation  (left)  to  a  4-class 
segmentation  (right)  of  a  rural  scene. 


Another  approach  to  layered  segmentation  is  to  first  segment  the  image 
into  regions  representing  various  textures  and  a  further  category  representing 
constant  tones  (no  texture).  Within  the  nontextured  regions  we  perform  a  seg¬ 
mentation  based  on  gray  level.  The  results  are  then  combined  to  produce  the 
final  segmentation  of  the  image.  This  technique  is  particularly  useful  for  the 
following  reason.  In  segmenting  an  image  for  texture  it  is  desirable  to  remove 

large  overall  variations  in  tonality  and  render  the  image  an  overall  middle 

1  2 

gray.  *  This  eliminates  artifacts  in  the  segmentation  due  to  variations  in 
illumination,  shadows,  and  so  forth.  However,  preprocessing  in  this  manner 
removes  tonal  features  in  nontextured  regions  that  may  be  important  for  a 
complete  and  meaningful  segmentation  of  the  image.  A  meaningful  result  can 
be  obtained  by  using  a  layered  approach  where  a  texture -oriented  segmentation 
of  the  preprocessed  image  is  followed  by  a  gray-tone  segmentation  of  the  un¬ 
processed  image. 

Figure  2  illustrates  this  approach  for  a  scene  containing  trees,  water, 
grass  or  fields,  and  a  feature  (apparently  a  road  or  bridge)  passing  over  a 
small  stream.  Figure  2  (upper  left)  shows  the  original  image  and  Fig.  2  (upper 
right)  shows  the  preprocessed  version.  Note  that  in  the  proprocessed  image 
the  bright  white  area  of  the  bridge  has  approximately  the  same  tonal  value  as 
the  dark  shadow  under  the  bridge.  Figure  2  (lower  left)  allows  the  result  of 
segmenting  the  original  image  into  regions  with  three  tonal  values  (dark,  shown 
in  red,  middle  grays,  shown  in  green,  and  light,  shown  in  blue).  Figure  2 
(lower  right)  shows  the  result  of  segmenting  the  preprocessed  image  into  tex¬ 
tures  representing  the  trees  (yellow),  the  grass  and  field  areas  (purple),  and 
the  nontextured  regions  such  as  water,  the  bridge,  and  the  shadows  (light  blue). 
Figure  3  shows  the  original  image  with  segmentation  boundaries  overlaid  (be¬ 
tween  tree  and  field  regions  and  between  trees  and  water)  in  solid  lines.  In 
addition,  the  boundary  obtained  by  overlaying  the  tonal  separation  of  Fig.  2 
(lower  left)  on  the  nontextured  areas  of  Fig.  2  (lower  right)  is  shown  in  gray. 
This  permits  us  to  obtain  a  more  accurate  representation  of  the  bridge  segment. 

As  another  example  of  this  layered  approach,  consider  the  portion  of  the 
scene  shown  in  the  white  box  in  Fig.  4(a)  and  enlarged  in  Fig.  4{b).  The  en¬ 
larged  scene  is  256  x  256  pixels  in  slice  and  represents  the  full  resolution  of 
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Fig.  2.  Layered  approach  to  segmenting  a  rural  scene  with  textured 
and  nontextured  areas.  Original  image  (upper  left),  preprocesaed 
image  (upper  right),  texture  segmentation  (lower  left),  and  gray-* 
level  segmentation  (lower  right). 
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Fig.  3.  Rural  scene  with  boundaries 
resulting  from  layered  segmentation 
superimposed. 


Fig.  4.  Enlarged  portion  of  a  river  scene  used  for  layered  segmentation 
experiment. 


the  data.  (The  photograph  from  which  the  digital  data  were  taken  was  scanned 
at  a  resolution  of  40  pm.)  Figure  5  shows  the  original  image  and  the  result 
aftei  segmenting  the  image  into  textured  regions  (trees)  and  nontextured  re¬ 
gions  and  applying  three-level  gray-tone  segmentation  within  the  nontextured 
regions.  Note  that  the  dark  shadows  that  the  trees  cast  on  the  water  are  easily 
separated  from  both  the  trees  and  the  water  using  this  procedure.  The  isola¬ 
tion  of  shadowed  areas  may  be  important  in  deducing  the  height  of  objects  in  a 
scene. 

A  final  example  of  segmenting  textured  from  nontextured  regions  is  shown 
in  Fig.  6.  The  original  scene  of  roadways  connecting  missile  sites  is  shown  in 
the  (a)  part  of  the  figure.  The  two-category  segmented  result  is  shown  in  the 
(b)  portion.  The  segmentation  of  this  scene  is  a  difficult  problem  because  the 
roads  tend  to  bleed  into  the  surrounding  terrain  in  many  places  and  the  pres¬ 
ence  of  features  on  the  roads  (3uch  as  the  objects  in  the  corner  of  the  parking 
lot)  tend  to  produce  something  like  a  small  textured  area.  Nevertheless,  a 
fairly  good  approximation  to  the  roadways  is  obtained.  Ideally  we  would  wish 
to  follow  this  texture  segmentation  with  a  gray-level  segmentation  to  highlight 
features  in  the  road.  However,  these  features,  because  of  their  size,  have 
tended  to  be  confused  with  the  texture  and  thus  many  do  not  appear  in  the  re¬ 
gion  designated  by  the  algorithm  as  roadways.  This  is  a  problem  that  requires 
additional  work.  The  features  in  the  roadway  can  perhaps  be  detected  by 
scrutinizing  the  error  residuals  in  the  segmentation  algorithms  or  possibly  by 
employing  more  texture  classes. 

Our  current  research  is  directed  toward  resolving  some  of  the  problems 
described  above  and  to  developing  a  method  for  automatic  selection  of  the  train¬ 
ing  data  (by  the  segmentation  algorithm).  The  latter  would  allow  the  segmenta¬ 
tion  algorithm  to  function  in  a  nonsupervised  mode  and  eliminate  the  need  for 
human  interaction  to  define  sets  of  training  data. 

3.  TECHNIQUES  FOR  REMOVING  DEGRADATIONS  CAUSED 

BY  LIGHT  CLOUD  COVER 

In  this  section,  we  will  discuss  the  use  of  homomorphic  filtering  to  expose 
objects  beneath  light  cloud  cover.  In  particular,  a  deterministic  model  of 
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Fig.  5.  Segmentation  of  the  river  scene  in  Fig.  4(b). 
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imaging  above  cloud  cover  motivates  an  approach  which  utilizes  an  adaptive 
homomorphic  filter.  This  space -varying  filter  is  parameterized  by  the  local 
mean  level  which  reflects  the  degree  of  the  local  degradation. 

The  approach  we  take  is  distinctly  different  from  other  homomorphic  filter¬ 
ing  procedures  for  image  enhancement  or  restoration.  It  differs  from  the  work 

3 

of  Mitchell  and  Delp,  which  recovers  images  degraded  by  light  cloud  cover, 

since  their  scheme  relies  on  a  stochastic  image.  It  differs  from  the  deter- 

4  5 

ministic  approach  of  Oppenheim  et  aL  and  Gilkes  for  image  enhancement  of 

cloudless  images,  since  again  theirs  are  nonadaptive  procedures.  In  essence, 

6  5 

it  is  closer  to  the  adaptive  approaches  of  Peli  and  Lim  and  Gilkes  where  an 
adaptive  filter  is  parameterized  by  the  local  deterministic  characteristics  of 
the  data. 

One  long-space  model  of  a  cloudy  image  is  described  in  Ref.  3.  Specifi¬ 
cally,  this  model  of  a  cloudy  image  is  stochastic  and  is  a  product  of  the  cloud 
transmission  function  and  a  function  of  the  ideal  image.  Our  cloudy  image 
model  is  deterministic  and  applies  on  a  short-space  basis.  It  assumes  that 
the  logarithm  of  the  image  can  be  divided  into  two  approximately  disjoint  spec¬ 
tral  bands:  the  cloud  spectrum  occupying  low  frequencies  and  the  image  spec¬ 
trum  occupying  high  frequencies.  Although  the  image  contains  some  low- 
frequency  information,  we  shall  assume  it  is  not  significant  in  exposing  object 
shapes  under  light  cloud  cover.  In  addition,  we  depart  from  a  typical  assump¬ 
tion  that  the  desired  image  itself  can  be  modeled  as  the  product  of  illumination, 
a  low-frequency  component,  and  reflectivity,  a  high-frequency  component. 
Rather,  we  view  the  illumination  component  in  the  same  way  as  we  view  the 
reflectance  component,  i.e.,  as  having  an  important  high-frequency  component 

7 

due  to  the  interaction  of  light  and  ground  objects.  Furthermore,  we  assume 
this  high-frequency  component  is  approximately  disjoint  from  the  cloud  transfer 
function  which  has  a  low -pass  characteristic. 

Before  proceeding  with  the  development  of  our  new  methods  and  compari¬ 
sons,  we  review  some  Important  ideas  and  formulate  a  framework  for  our 
investigations. 
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3.1  Modeling  the  Imaging  Process 

A  number  of  approaches  to  modeling  the  imaging  process  have  been  pre¬ 
sented  in  the  literature.  A  more  difficult  problem  is  to  model  the  imaging 
process  with  light  cloud  cover.  In  this  section,  we  first  review  one  approach 
to  modeling  undegraded  images.  Two  different  viewpoints  are  presented  which 
rely  on  an  illumination-reflectivity  model.  We  then  enter  a  stochastic  frame¬ 
work  in  which  a  model  of  images  degraded  by  light  cloud  cover  will  be  dis¬ 
cussed.  Finally,  we  present  our  own  deterministic  interpretation  of  cloudy 
images. 

In  the  formation  of  images,  the  illumination  and  reflectivity  are  combined 
4 

by  a  multiplication  law: 

I(n,  m)  =  i(n,  m)  •  r(n,  m)  (2) 

where  i(n,  m)  and  r(n,  m)  are  the  illuminance  and  reflectance  components,  re¬ 
spectively.  One  assumption  is  that  the  illumination  varies  slowly  (i.e.,  it  con¬ 
tains  mainly  low-frequency  components)  and  the  reflectance  is  sometimes  dy¬ 
namic  and  sometimes  static  and  therefore  may  be  regarded  as  containing  mainly 
high-frequency  components.  The  logarithm  operation  separates  the  multiplica¬ 
tive  signal  into  two  additive  components : 

log  (I(n,  m)J  o  log (i(n,  m)]  +  log  (r(n,  m)]  .  (3) 

If  a  linear  amplifier  or  attenuator  with  gain  a  follows  the  log,  the  image  out¬ 
put  is  given  by 

l'(n,m)  s  {I(n.m))0'  .  (4) 

When  a  <  i,  we  get  a  washed-out  image  and  when  a  >  1,  the  image  is  sharpened 
but  might  be  saturated.  If  we  want  to  reduce  the  dynamic  range,  which  Is  con¬ 
tributed  mostly  by  the  illumination,  and  increase  the  sharpness,  which  is  con¬ 
tributed  mostly  by  the  reflectance,  we  should  choose  a  as  a  function  of  fre¬ 
quency,  i.e.,  a  <  i  for  high  frequencies  as  in  Fig.  7.  This  is  the  essence  of  the 
homomorphic  filtering  operation. 

Since  the  illumination  is  not  exclusively  low  frequency  and  the  reflectance 
consists  of  both  high-  and  low-frequency  components,  the  output  of  the  homo¬ 
morphic  filtering  operation  has  some  artifacts.  It  seems  that  the  degree  to 


20 


which  those  artifacts  are  visible  is  strongly  dependent  on  the  way  a  changes 
from  a  <  1  to  a  >  1.  As  pointed  out  by  Schreiber/  homomorphic  filtering 
can  be  replaced  by  a  linear  filter  followed  by  a  power  law  for  low-contrast 
images.  Schreiber  hypothesizes  that  in  natural  scenes,  the  illumination  func¬ 
tion  is  as  essential  to  perception  as  the  reflectivity  of  the  object.  Both  have 
a  broad  spectrum,  generally  with  more  power  at  the  low  spatial  frequencies, 
The  reflectance  derives  its  low-frequency  component  from  the  presence  of 
large  patches  of  relatively  constant  value.  The  illumination  derives  its  high- 
frequency  component  from  interaction  between  the  edges  and  surfaces  of  ob¬ 
jects,  which  are  at  many  different  angles,  and  the  incident  light.  Consequently, 
any  attempt  to  separate  these  components  by  homomorphic  filtering  will  be 
limited  in  its  success,  instead,  adaptive  control  of  the  low  and  high  components 
of  the  original  signal  is  preferable  for  enhancement. 

Now  let  us  examine  one  model  of  images  degraded  by  cloud  cover.  In  this 
model  the  ground  reflection  passes  through  the  clouds.  The  energy  collected 
above  the  cloud  consists  of  two  components :  ground  information  and  scatter¬ 
ing  from  the  clouds.  Figure  8  is  a  simplified  model  of  the  Imaging  process. 
The  energy  recorded  is  given  by: 

s(n,  m)  =  aLr(n,  m)  t(n,  m)  +  L  [1  -t(n,  m)]  (5) 

where 

L  =  sun  illumination 

a  =  attenuation  of  the  illumination  in  a  downward 
direction  (assumed  constant) 
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r(n,  m)  =  reflectance  of  the  ground 

t(n,m)  =  transmittance  of  the  cloud  (upward), 

0  <  t(n,  m)  <  1 

In  order  to  recover  the  desired  reflectance  information,  it  was  assumed 
that  the  transfer  function  of  the  cloud  and  the  reflectance  of  the  ground  are 
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Fig.  8.  A  simplified  model 
of  imaging  through  clouds. 


stochastic  processes.  It  was  also  assumed  that  the  transmittance  of  the 
cloud  has  relatively  more  energy  in  low  spatial  frequencies.  From  Eq.  (5), 
we  can  obtain  a  multiplicative  form  of  the  '’noise"  and  the  "signal"  given  by: 

L  -  8{n,  m)  =  t(n,  m)  (L  -  aLr(n, m}]  .  (6) 

Hy  taking  the  log  of  Eq.  (6),  we  obtain  tire  sum 

log(L  -  s(n,  ml!  =  log(t(n,  m)l  +  log{ L  -  aLr{n,  m)l  (7) 

or 

P(n,  m)  =  N(n,m)  +  M(n,  m) 

where 

P(n,  in)  =  log  (L  -  s(n,  m)]  =  "signal  +  noise" 

N(n,m)  =  log(t(n,  m)J  *  "noise" 

M(n,  m)  =  log[L  -  aLr(n,m))  =  "signal" 
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3.2 


Filtering  Based  on  the  Image  Model 


3 


The  recovery  of  the  signal  from  Eq.  (7)  can  be  approached  by  the  applica¬ 
tion  of  Wiener  filtering.  If  we  assume  that  the  signal  and  the  noise  are  uncor¬ 
related,  we  can  reduce  the  noise  by  filtering  the  given  degraded  image  with 
the  following  optimal  filter 


. .  .  +  V  Mp-KW 

'  2’  =  +  SNN<‘V“2>  +  2V 


'  yv? 


(8) 


where 


M  is  the  mean  value  of  the  given  signal  M 
is  the  mean  value  of  the  given  noise  N. 


In  order  to  estimate  Spp,  the  power  spectrum  of  the  signal  plus  noise  and 
Snn.  the  power  spectrum  of  the  noise,  the  values  L  and  t(n,  m)  need  to  be 
calculated.  L  was  estimated  as  the  highest  value  in  the  image  and  t(n,  m)  was 
estimated  by  the  following: 


A 

t(n,  m) 


L  -  S(n,  m) 

~r-e— 


(9) 


where  G  is  a  constant  which  approximates  the  typical  ground  reflection 
(i.e.,  the  average  of  aLr(n,  m)]. 

The  power  spectrum  of  the  noise  was  estimated  roughly  by  the  magnitude 
squared  of  the  Fourier  transform  of  t(n,  m).  Since  the  transfer  function  of  the 
cloud  was  considered  to  possess  only  low  frequencies,  the  estimate  of  this 
power  spectrum  was  low -pass  filtered.  The  cross  section  of  the  resulting 
filter  is  given  in  Fig.  9. 

The  filter  has  a  boost  at  (w^,  =  (0,  0)  as  a  result  of  assuming  a  non¬ 

zero  mean  process  and  an  attenuation  at  the  low  frequencies.  Except  for  the 
value  at  (0,  0),  the  filter  has  basically  the  same  functional  form  as  the  gain  in 
the  homomorphic  filtering  approach.  It  should  be  noticed  that  in  homomorphic 
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|  H  (u, ,  0)  I 


Fig.  9.  Cross  section  of  the  2-D 
Wiener  filter. 


filtei*ing,  it  is  possible  to  sharpen  the  edges  while  in  the  Wiener  filter  approach, 
it  is  not  since  |H(u>)|  <>  1  for  all  w  ^  0. 

In  general,  the  noise  spectrum  density  is  not  known  accurately 

and  consequently,  the  Wiener  filter  used  to  process  the  noisy  image  may  be 
suboptimal.  An  alternative  approach  involves  iteratively  updating  the  esti¬ 
mate  of  the  spectral  density  of  the  noise.  Since  an  estimate  of  the  reflectance 
function  aLr(n,m)  is  used  in  determining  (3  in  Eq.  (9),  the  improved  estimate 
of  aLr(n,  m)  obtained  from  the  Wiener  filter  output  can  be  used  to  update  the 
estimate  of  the  noise  spectral  density.  The  Wiener  filter  in  Eq,  (8)  is  then  re¬ 
computed  and  the  original  cloudy  image  is  again  filtered.  If  we  do  this  re¬ 
peatedly,  the  resulting  iterative  process  should  converge  to  a  better  estimate 
than  generated  by  the  one-step  process.  In  the  space  domain,  the  iterative 
cloud  estimate  is  given  by 


(n,  m) 


L  -  s(n,  m) 

-TO— 


(40a) 


where  Gq  is  a  constant.  From  ^(n,  m),  we  get  a  better  estimate  of  aLr(n,  m) 
denoted  by  G^in,  m),  and  the  iteration  is  continued  as 


(n,  m) 


L  -  S(n,  m) 
L-G, 


where  Gj. 


is  the  kth 


estimate  of  the  function  aLr(n,  m). 


(10b) 
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3.3  Homomorphic  Filtering 

The  assumption  that  images  and  clouds  are  stochastic  processes  is  very 

artificial,  A  more  appropriate  approach  invokes  the  assumption  that  the  two 

7 

are  deterministic  processes.  In  fact,  we  adopt  Schrelber's  approach  which 
assumes  that  both  illumination  and  reflectance  contain  low  and  high  components 
and  that  the  high  frequencies  are  more  Important  for  perception.  At  the  same 
time,  we  assume  that  the  transfer  function  of  the  cloud  contains  only  low 
spatial  frequencies.  We  now  interpret  Eq. (6)  so  that  tin,  ir)  is  determin¬ 
istic  and  contains  mainly  low  frequencies  and  L  -  aLrin,  m)  is  also  determin¬ 
istic  and  contains  perceptually  important  high-frequency  information.  8y 
using  these  assumptions,  the  image  aLr(n,  m)  can  be  restored  by  homomorphic 
filtering. 

As  discussed  in  Sec.  3,2,  both  the  image  illumination  and  reflectance  con¬ 
tain  low-  and  high-frequency  components.  In  our  application,  however,  the 
perceptually  important  high-frequency  component  of  illumination  is  actually  a 
blessing.  That  is,  since  light  cloud  cover  is  primarily  low  pass,  it  can  be 
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homomorphically  filtered  without  significantly  disturbing  the  high-frequency 
component  of  the  desired  image's  illumination  and  reflectance.  Of  course,  the 
desired  low-frequency  components  will  be  altered.  Nevertheless,  these  changes 
should  reflect  themselves  simply  as  illumination  changes  over  large  constant  re¬ 
gions,  thus  not  disturbing  any  significant  information  under  the  light  cloud  cover. 

An  implied  assumption  here  is  that  the' cloud  spectrum  contains  predomi¬ 
nantly  low-frequency  components,  which  motivates  the  long-space  homomor¬ 
phic  filtering  solution.  In  fact,  this  is  similar  to  the  assumptions  made  by 
3 

Mitchell  et  al.  in  a  stochastic  framework.  However,  the  thickness  of  cloud 
cover  can  change  over  the  extent  of  an  image.  Consequently,  an  adaptive 
approach  should  be  better  suited  to  restoring  the  desired  image,  i.e.,  an  ap¬ 
proach  where  the  homomorphic  filter  changes  with  local  cloud  characteristics. 

6  8 

In  earlier  work,  *  we  developed  a  locally  adaptive  contrast  enhancement 
technique  based  on  high-pass  filtering.  The  local  mean  of  the  cloudy  image 
was  used  to  determine  how  much  gain  was  applied  to  the  high-frequency- com¬ 
ponents  of  the  image  in  a  particular  area.  In  this  way  the  amount  of  contrast 
increase  was  adaptively  adjusted  across  the  entire  image.  The  adaptive  homo¬ 
morphic  filtering  approach  described  here  is  similar  in  that  measurements 
made  on  local  image  characteristics  are  used  to  determine  the  homomorphic 
filter  characteristics. 

The  adaptive  technique  can  be  formulated  as  follows.  Consider  applying 
a  short-space  window  to  the  noisy  signal  Pin,  m)  in  Bq.  0): 

Pj  ^(n,  m)  *  Wf  j.(n,  m) »  Pin.  m)  -  W(n  -  I  M/2,  m  -  kM/2) 

•  P(n,m)  .  (12) 

Let  us  assume  the  window  takes  on  a  pyramidal  shape  and  is  shifted  over  the 
data  at  intervals  of  half  its  length  (i.e.,  M/2  where  M  x  M  is  the  extent  of  the 
window).  Consequently,  the  window  1ms  the  following  desirable  property: 

£  £  Win  -  I  M/2,  m  -  kM/2)  =  1  .  (13) 

/  k 

Applying  an  adaptive  filter  which  operates  on  each  segment  j.(n.  m),  we 
obtain: 
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Pi)k(n,m)  =  P^k(n,m)  *+  h^k(n,m,0) 

-  [W(n-iM/2,  m-kM/2)  P(n,m)]  **  h|  k(n,m,0)  (14) 

where  the  parameterized  flPer  impulse  response  h^  k(n,m,0)  is  a  function  of 
the  vector  0,  which  is  a  set  of  parameters  dependent  on  the  local  cloud  char¬ 
acteristics.  (The  double  asterisk  denotes  2-D  spatial  convolution.)  In  partic¬ 
ular,  0  is  a  function  of  the  DC  level  of  the  windowed  signal  which  reflects  the 
cloud  density  under  the  window.  Further,  the  filter  is  a  high  pass  where  the 
shape  and  amplitude  depend  on  0 .  Later  we  will  elaborate  on  this  design. 

For  Eq.  (7),  the  filtering  process  is  then  given  by 

Pf  k(n,  m)  =  |w(n  -  i M/2,  m  -  kM/2)  {log  (t(n,  m)] 

•f  log{L  -  aLr(n,  m»)  j  **  h^  k(n,m,0)  .  (15) 

We  assume  that  the  window  is  "sufficiently  smooth"  so  that,  the  low-pass 
nature  of  the  noise  term  log  (t(n,m)|  and  the  high-pass  nature  of  the  sJgn«l 
component  lc^  (L  -  aLr(n,  m)]  are  preserved  after  windowing.  Consider  the 
case  wh-,  re  the  noise  and  the  signal  are  disjoint  in  frequency,  and  the  filter 
itf  K(n,  m,0)  is  an  ideal  high-pass  filter  whose  nonzero  energy  band  matches 
that  of  tlx©  signal,  then  we  can  write 

P«  .  (n,  m)  a  W(n  -  iM/2,  m  -  kM/2)  log  (L  -  aLr(n,  m))  .  (16) 

In  practice,  however,  these  assumptions  do  not  hold  exactly.  Since  the  noiao 
and  signal  are  only  approximately  disjoint,  the  noise  term  will  not  be  entirely 
removed,  and  since  the  filter  is  not  an  ideal  high  pass,  the  signal  will  not  re¬ 
main  intact.  Alternatively,  whe.  k(n,  m,  0)  adaptively  amplifies  the  high 
frequencies  and  attenuates  the  lows,  the  cloud  will  be  only  partially  suppressed 
{while  somewhat  disturbing  the  low-froquency  component  of  log[L  -  al,  r(n,  m)]} 
and  the  high-frequency  detail  of  the  desired  signal  will  be  enhanced  (while 
boosting  any  cloud  energy  in  this  region).  The  relation  in  Eq-  (16)  represents 
an  approximation  of  the  desired  windowed  signal  along  with  the  effects  of  any 
residual  cloud  noise.  The  parameters  0,  in  hf  k(n,  m,0),  which  rely  on  the 
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local  cloud  density  should  be  chosen  so  that  the  noise  is  suppressed  as  much 
as  possible  without  altering  the  desired  signal. 

With  these  approximations  in  mind,  our  reconstruction  procedure  (which 
we  refer  to  as  ovcrlap-add)  can  then,  with  the  use  of  Eq.  (13),  be  written  as 

P(n,  m)  =  /]  Yi  W(n  ~  Mk/2,  m  -  Mk/2)  log[L  -  aLr(n,m)] 
i  k 

=  log[L  -  aLr(n,  m)]  £  J]  W(n  -  Mi/2,  m  -  Mk/2) 

i  k 

=  log  [L  -  aLr(n,  m)]  .  (i7) 

This  procedure  is  illustrated  in  Fig.  10.  Finally,  to  obtain  an  estimate  of  the 
signal,  we  exponentiate  P(n,m)  to  recover  aLr(n.m). 


Fig.  10.  Overlap-add  technique  for  adaptive  image  filtering. 
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From  the  simple  model  of  imaging  above  clouds,  we  have  concluded  that 
it  is  desirable  to  adaptively  attenuate  the  low  frequencies  and  to  amplify  the 
high  frequencies  of  an  image.  A  filter  which  has  the  desirable  shape  and  is 
computationally  efficient  is  a  circularly  symmetric  Gaussian  which  is  shifted 
to  w  =  ir.  A  cross  section  of  this  filter  is  given  in  Fig.  11  where  |  H(0)  |  <  1, 


H(w1,o)2)  =  A 


(18) 


The  filter  size  was  chosen  to  be  fixed  at  16  x  16.  The  filter  parameters  A,  B, 
and  C  are  functions  of  the  value  of  the  Fourier  transform  of  the  windowed  data 
at  (w.,w2)  =  (0,0)  (i.e.,  the  DC  value  of  the  pyramid -shaped  window).  The 
filter  parameters  are  given  for  two  extreme  DC  values:  (Al,  Bl,  Cl)  for  a 
DC  level  of  zero  and  (A2,  B2,  C2)  for  a  DC  level  of  225.  For  a  DC  value  that 
lies  between  these  two  extremes,  the  parameters  are  computed  by  a  quadratic 
equation 

2 

(Y  -  Y  )  •  D 

Y  =  — - - K -  +  Y.  for  Y  =  A,  B,  C  (19) 

255*  1 

where  D  denotes  the  DC  level  and  where  Y^  and  Y2  represent  the  two  extreme 
values  of  each  parameter. 

It  is  desirable  that  for  a  low  value  of  the  DC  level  (i.e.,  little  or  no  cloud 
cover),  very  little  attenuation  and  amplification  should  be  used  and,  further, 
the  amplification  should  be  applied  to  only  the  very  high  frequencies  (i.e.,  B, 
the  standard  deviation  of  the  Gaussian  shape  in  Fig.  11,  should  be  small).  For 
a  high  DC  level  (i.e.,  thick  clouds),  the  attenuation  and  amplification  should 
be  applied  to  a  wider  range  of  high  frequencies  (i.e.,  B  should  be  larger).  A 
quadratic  interpolation  for  computing  A,  B,  and  C  was  used  rather  than  a 
linear  interpolation  because  a  larger  range  of  low  level  of  luminance  needs 
almost  no  processing. 
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Fig.  11.  Cross  section  of  a  Gaussian 
high-pass  filter. 


3.4  Experimental  Results 

An  experiment  was  performed  in  order  to  compare  the  Wiener  filtering 
approach  to  the  homomorphic  filtering  approach  for  removing  degradations 
caused  by  light  cloud  cover.  An  aerial  reconnaissance  image  partially  de¬ 
graded  by  light  cloud  cover  is  shown  in  Fig.  12.  Figure  13  is  the  result  of  the 
Mitchell  and  Delp1  algorithm  when  the  estimate  of  the  typical  ground  reflection 
G  was  chosen  to  be  a  constant  equal  to  140.  The  processed  image  seems  to 
contain  only  the  high  frequencies  of  the  unprocessed  image. 

By  using  the  iterative  approach  described  by  Eqs.  (10a)  and  ( 1  Ob)  one  can 
improve  the  estimate  of  the  noise  spectrum.  Figure  14  contains  the  result 
after  three  iterations. 

Figure  15  is  the  result  after  three  iterations  of  improving  the  estimate 
of  the  noise  and  also  reprocessing  the  processed  image  as  in  Eq.  (ila).  One 
can  see  that  the  second  form  of  iteration  is  more  desirable. 

Long-space  homomorphic  filtering  was  examined  with  two  different  filters. 
Each  lias  a  Gaussian  shape  as  described  in  Sec.  3.3  and  a  length  of  256  x  256 
(i.e,,  the  length  of  the  unprocessed  image).  The  parameters  of  the  first  filter 
are  A  =  3,65,  B  =  320,  and  C  =  -2.15.  The  processed  image  is  shown  in 
Fig.  16.  The  second  filter  which  sharpens  and  attenuates  more  than  the  first 
lias  the  parameters  A  =  6,21,  B  =  320,  and  C  *  -4.21.  The  processed  image 
is  illustrated  in  Fig.  17.  From  these  results,  we  can  see  that  there  is  some 
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Aerial  photograph  with  some  light  cloud  cover. 


honoi-j! 

•  '  *>  - 

" 

'  '  ' 

•  'v. • 

,'V\  , 
v  .\x  ^  ■ 

;  • i  i% 

■a 

.  C  '  i 

■  MV'  ■ 

w  * 

Fig.  13.  Image  processed  by  application  of  a  global  Wiener  filter 
applied  to  log(L  -  s(n,m)J.  (See  Ref.  1.) 


1 109703- s] 


Fig.  14.  Image  processed  with  a  global  Wiener  filter  iterated 
three  times  to  obtain  improved  estimate  of  noise  power 
spectrum. 


Fig.  15.  Image  processed  according  to  Eq.  (11)  with  three  iterations. 


Fig.  16.  image  processed  with  a  long-space  homomorphic  high-pass  filter 
with  parameters  A  .  3.65,  B  =  320,  and  C  .  -2.15. 
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Fig.  17.  Image  processed  with  a  long-space  homomorphic  high-pass  filter 
with  parameters  A  9  6,21,  B  =  320,  and  C  =  -4.21. 


sharpening  but  most  of  the  cloud  was  not  removed  and,  therefore,  some  ob¬ 
jects  were  not  exposed. 

In  the  adaptive  homomorphic  filtering  approach,  a  pyramidal  window  of 
size  16  x  16  was  chosen  and  the  parameters  ranges  are  given  in  Table  I.  The 
filter  shapes  are  shown  in  Fig.  18,  and  the  processed  image  is  illustrated  in 
Fig. 19. 


TABLE  1 

RANGE  OF  PARAMETERS  FOR  THE  ADAPTIVE 
HOMOMORPHIC  FILTERING 

DC 

0 

255 

A 

0.2 

1.82 

B 

5.0 

20.0 

C 

0.9 

-0.32 

Fig.  18.  Range  of  filter  shapes  for  adaptive  homomorphic  system. 


We  see  from  our  results  that  the  most  desirable  restoration,  i.e.,  maxi¬ 
mum  suppression  of  the  cloud  and  enhancement  of  the  objects,  was  achieved 
by  adaptive  homomorphic  filtering.  We  are  beginning  to  explore  implementa¬ 
tion  strategies  for  this  algorithm  to  improve  its  computational  efficiency. 

4.  ERRORS  CAUSED  BY  TRUNCATION  EFFECTS 
IN  THE  ITERATIVE  IMPLEMENTATION  OF 
MULTI-DIMENSIONAL  DIGITAL  FILTERS 

When  an  iteration  is  used  to  compute  successively  better  approximations 
to  the  output  of  a  digital  filter,  the  effective  impulse  response  grows  in  extent 
at  each  iteration.  For  the  case  where  the  frame  storage  buffer  has  a  fixed 
size,  tiie  extent  of  the  successive  output  estimates  is  limited  by  the  buffer  size, 
resulting  in  spatial  truncation  errors  which  affect  successive  approximations 
to  the  desired  output  signal.  In  this  section,  we  shall  discuss  our  preliminary 
findings  on  truncation  errors  and  their  relation  to  boundary  conditions.  We 
begin  by  giving  a  brief  review  of  the  iterative  implementation. 


4.1  Brief  Review  of  the  Iterative  Implementation 

The  iterative  filter  possesses  an  interesting  structure  to  implement,  since 
it  involves  FIR  (finite -extent  impulse  response)  filtering  operations  as  well  as 
the  feedback  of  buffered  output  image  frames.  Because  of  this  frame  feedback, 
the  iterative  filter  can  be  used  to  approximate  the  output  of  a  2-D  symmetric 
rational  transfer  function.  Consider  the  filter  frequency  response  of  the  form 

A(w4,  w-) 

(20) 

where  A(Wj,u>2)  and  PfW  j.Wg)  are  2-D  trigonometric  polynomials  with  coeffi¬ 
cients  altij.iij)  and  respectively.  With  appropriate  normalization, 

we  can  write 

B(w1,u>2)  *  i  -  C(a.'1,o)2)  (21) 

where  C(W|,w2)  Is  another  trigonometric  polynomial  with  coefficients  c(n1(n2). 
Let  ua  filter  an  input  image  x(n1(  n2)  to  obtain  an  output  Image  y(n1#  «2).  Then 
x{nlt  n2)  and  y(nt,n2)  satisfy  the  implicit  relation 
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(22) 


y(n1#n2)  =  a(nltn2)  **  x(nd,n2)  +  c(n1#n2)  **  y(n1,n2) 

where  the  double  asterisk  denotes  2-D  convolution.  We  can  use  this  relation 
iteratively  to  generate  successively  better  approximations  to  y(n^,  n2 ).  Letting 
i  denote  the  iteration  number,  we  can  write 

yi<nl,n2)  =  a(n14n2)  **  xi(n1>n2)  +  c(n1,n2)  **  (23a) 

or  in  the  frequency  domain 

Yil«i,«2)  =  A(uitu>2)  +  C(w1,«2)  Yi_i(a»i,w2)  .  (23b) 

The  subscript  i  on  the  input  image  x.(n^,n2)  is  included  to  indicate  that  the 
iterative  implementation  may  be  generalized  to  accept  sequences  of  image 
frames  as  its  input.  Similarly,  the  image  sequence  y,(nj,n2)  may  be  regarded 
as  a  succession  of  image  frames.  At  this  point,  it  is  natural  to  consider  the 
iteration  index  i  as  a  time  variable. 

Figure  20  shows  a  simple  block  diagram  for  the  iterative  implementation 
which  highlights  the  basic  operations.  There  are  two  FIR  filtering  operations, 
denoted  by  the  boxes  labeled  A(e>j,c»>2)  and  CfWj.Wg),  a  summing  node,  and  a 
frame  storage  buffer. 


Fig.  20.  Block  diagram  of  the  2-D  iterative  filter  implementation. 
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4.2  Characterization  of  the  Truncation  Error 


The  operation  of  truncation  is  particularly  important  since  we  cannot  deal 
in  practice  with  infinitely  long  sequences.  Because  the  output  of  a  rational 
filter  has  infinite  extent  in  general,  errors  will  be  introduced  by  repeated 
truncation  and  may  become  intolerably  large  over  the  limited  domain  of  the 
output  signal.  The  error  over  this  limited  region  can  be  viewed  as  the  solu¬ 
tion  to  a  homogeneous  partial  difference  equation  with  specified  boundary  con¬ 
ditions.  In  particular,  the  error  is  linearly  dependent  on  values  of  the  ideal 
solution  along  the  boundary  of  the  region  remaining  after  truncation. 

The  iteration  (23a)  can  be  written  in  operator  notation  as 

yi(n1,n2)  =  Fy._1(n1,n2)  (24) 

where  F  is  an  operator  of  the  form 

Fy(n1#n2)  =  a(n1,n2)  x(n1,n2)  +  ctr^,^)  **  y(nt,n2)  .  (25) 

Under  the  appropriate  assumptions,  there  exists  a  unique  solution  y(n^,  n2)  = 
Fy(n^,  n2).  When  the  size  of  the  frame  buffer  is  exceeded,  the  2-D  signal 
yi(nl  *  n2 )  must  ^e  truncated  in  extent.  This  introduces  a  truncation  operator 
T  into  the  iteration,  giving  Eq.  (24)  the  form 

yi(n1,n2)  =  TFyl_1(n1,n2)  (26) 


where 


Ty(n1,n2)  = 


y(n1#  n2) 


0 


for  (nltn2)  gl 
for  (nlfn2)  f*I 


(27) 


(The  region  I  can  ba  thought  of  as  representing  the  extent  of  the  frame  buffer.) 
The  iteration  (26)  also  can  be  shown  to  have  a  unique  solution  y(n,,n_)  * 

A  1 

TFy(n1,n2).  In  general,  the  solution  to  Eq.  (24),  yfn^n^,  will  not  be  equal 
to  y^.n-j).  The  difference  batween  these  two  signals  is  attributable  to  the 
truncation  effects.  Let  us  define  the  truncation  error  as 

e(n1.n2)=y(n1,n2)-y(n1,n2)  .  (28) 
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It  can  be  shown  that  the  error  signal  e(n^,n2)  satisfies  an  iteration  which 
corresponds  to  driving  Eq.  (23a)  with  no  input  [x^n^n^  =  0]  but  imposing 
boundary  conditions  around  the  edges  of  the  region  I.  Let  us  consider  a  re¬ 
gion  1^  which  surrounds  the  region  I.  The  extent  of  1^  depends  on  the  extent 
of  impulse  response  c(n1,  n2)  in  Eq.  (23a).  We  can  now  define  the  boundary 
conditions  bnd(n^,n2)  to  be 


bnd(n1#n2)  =  j 


yf^,^) 

o 


for  (n^,  n2)  €  Ib 
elsewhere 


(29) 


Then  the  error  signal  e(n^,n2)  satisfies  the  iteration 


ei(nr ~  TCc<ri1» n2>  **  ei-l(nl'n2^  +  ^d(n1#n2) 

for  (nrn2)  e  I  +  Ib  •  (30) 

The  energy  in  the  error  signal  can  be  shown  to  be  proportional  to  the  energy 
in  the  boundary  condition  signal  bnd(nj,  m,).  Thus,  if  the  sirie  of  the  frame 
buffer  (and  hence  the  extent  of  the  region  I)  is  large  enough,  one  would  expect 
that  the  correct  solution  y(nj,n2)  (and  hence  the  boundary  condition  signal 
bnd(n1(  »2)J  would  be  small  in  the  region  1^,  In  this  case,  the  error  signal 
efn^  n2)  also  will  be  small,  and  y(nj,  n2)  will  be  a  good  approximation  to  the 
desired  output  signal  y(n^,n2).  Furthermore,  the  error  can  be  eliminated  by 
including  the  boundary  condition  information  bnd(n^,n2)  in  the  spatially  trun¬ 
cated  iteration,  giving  it  the  form 

yl<nrn2>  3  TFyi-i(ni*n2)  +  bnd(ni*n2) 

for  (nrn2)  €  I  +  Ib  •  (31) 

5.  POTENTIAL  ARCHITECTURES  FOR  THE  ITERATIVE 
IMPLEMENTATION  OF  MULTI -DIMENSIONAL  DIGITAL 
FILTERS 

In  this  section,  we  shall  examine  some  architectures  which  support  the 
iterative  implementation  of  multi-dimensional  digital  filters.  We  begin  by 
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discussing  possible  data-handling  approaches  and  their  implications  with  re¬ 
spect  to  architectures.  Later,  we  shall  take  a  "software"  point  of  view  and 
discuss  the  problem  of  partitioning  signal-processing  operations  among  several 
processors  in  a  multiprocessor  architecture. 

Referring  back  to  Fig.  20,  we  see  that  the  iterative  implementation  con¬ 
sists  of  two  FIR  filtering  operations,  denoted  by  the  boxes  labeled  A(w^,w.>) 
and  a  summing  node,  and  a  frame  storage  buffer.  Consequently, 

much  of  the  discussion  of  processing  architectures  will  be  directed  toward 
the  efficient  implementation  of  multi-dimensional  FIR  filters. 

5.1  Image  Flow  Options 

There  are  three  popular  methods  for  describing  how  a  sequence  of  images, 
each  composed  of  many  picture  elements  (pixels),  may  flow  through  a  signal- 
processing  architecture.  The  iterative  implementation  as  diagrammed  in 
Fig.  20  suggests  a  "frame-by-frame"  image  flow  which  might  result  naturally 
from  a  lens  imaging  system;  all  the  image  pixels  from  a  given  frame  are 
available  simultaneously.  One  of  the  major  problems  with  this  image  flow  is 
providing  the  large  number  of  parallel  data  paths  necessary  to  shuttle  entire 
images  around  the  processing  structure.  Conversely,  one  advantage  is  the 
relatively  long  available  data  transfer  time,  typically  1/30  s.  Consequently, 
there  is  flexibility  in  trading  data  transfer  time  for  the  number  of  data  paths 
using  multiplexors  and  demultiplexors. 

At  the  other  extreme  of  the  image-data-flow  spectrum  we  have  the  "se¬ 
quential  pixel"  flow.  All  the  pixels  in  an  image  are  transferred  one  after 
another  (multiplexed),  followed  by  the  pixels  of  the  next  image  frame,  and  so 
on.  The  processing  structure  sees  a  continuing  stream  of  pixel  values.  This 
serial  flow  of  image  data  can  result  from  sensors  such  as  the  flying  spot 
scanner,  where  a  single  detector  is  swept  across  an  image  in  a  raster  scan. 
With  the  sequential  pixel  flow,  only  one  data  path  is  needed  since  all  the  pixels 
are  multiplexed  into  a  single  stream.  However,  the  data  transfer  rates  can 
become  quite  large.  For  a  high  resolution  1024  x  1024  image,  the  data  trans¬ 
fer  rate  is  approximately  30M  pixels/second. 


t 
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The  "parallel  row"  image  flow  offers  a  compromise  between  the  frame- 
by-frame  and  the  sequential  pixel  image  flows.  The  parallel  row  image  flow, 
as  the  name  implies,  transfers  rows  of  image  data  in  parallel,  but  the  pixels 
within  each  row  are  transferred  serially.  This  serves  to  reduce  the  number 
of  parallel  data  paths  compared  to  the  frame-by-frame  method,  and  yet  it  does 
not  require  the  high  transfer  rates  of  the  sequential  pixel  method.  Unfortu¬ 
nately,  the  number  of  parallel  data  paths  required  may  still  be  large  for  typical 
image  data.  The  parallel  row  image  flow  seems  potentially  well  suited  to 
sensor  systems  consisting  of  a  linear  array  of  detectors;  each  row  of  the 
image  corresponds  to  data  from  one  of  the  detectors. 

5.2  An  Architecture  for  Iterative  Filters  Using  Sequential 
Pixel  Image  Flow 

As  we  saw  in  Fig.  20,  the  iterative  filter  contains  two  2-D  FIR  convolu¬ 
tions,  so  let  us  begin  by  examining  the  structure  for  convolution  with  a  sequen¬ 
tial  pixel  image  flow.  The  approach  is  straightforward;  basically,  it  consists 
of  a  pixel  buffer  which  allows  simultaneous  access  to  all  the  input  pixels  neede< 
to  compute  a  particular  output  pixel.  Figure  21  shows  the  implementation  of 
a  3  x  3  FIR  filter  on  a  sequential  pixel  stream.  The  unit  delays  in  conjunction 
with  the  shift  register  delays  form  the  necessary  pixel  buffer.  The  length  of 
the  shift  register  is  two  less  than  the  length  of  an  image  row,  so  that  a  shift 
register  delay  plus  two  unit  delays  will  buffer  an  entire  row  of  pixels.  The 
appropriate  nine  pixels  are  multiplied  by  the  FIR  filter  coefficients  and  then 
summed  to  form  a  single  output  pixel.  The  gating  logic  shown  in  Fig.  21  allowE 
one  or  more  of  the  product  terms  to  be  zeroed  out;  this  is  necessary  to  handle 
edge  pixels  in  the  images  and  the  boundary  between  images. 

The  convolution  operation  shown  in  Fig.  21  may  now  be  used  in  the  imple¬ 
mentation  of  the  iterative  filter.  In  effect,  Fig.  21  is  Inserted  into  the  boxes 
labeled  and  C(c Fig.  20.  The  frame  storage  buffer  in  Fig.  20 

can  be  implemented  with  a  shift  register  whose  length  is  equal  to  the  number 
of  pixels  in  an  image.  It  may  even  be  possible  to  combine  the  two  filter  output 
accumulators  with  the  summing  node  in  Fig.  20  for  more  efficient  operation. 
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Fig.  22,  Structure  for  a  row-parallel  convolution.  (Braces  represent 
coefficient  multiplication,  accumulation,  and  gating.) 


There  are  some  flexibility  problems  with  this  approach.  The  size  of  the 
pixel  buffer  in  Fig.  21  depends  on  the  size  of  the  FIR  filter  being  implemented, 
as  does  the  number  of  coefficient  multipliers,  control  gates,  and  accumulator 
inputs.  A  structure  designed  to  implement  convolutions  of  a  fixed  size,  say 
3x3,  would  not  be  able  to  implement  larger-size  convolutions.  On  the  other 
hand,  a  processor  designed  to  support  an  11  x  11  convolution  would  not  be 
making  full  use  of  its  computational  potential  when  implementing  any  smaller- 
size  convolutions.  Furthermore,  the  sizes  of  tne  frame  storage  buffer  in 
Fig.  20  and  the  gating  logic  signals  in  Fig.  21  depend  on  the  number  of  pixels 
per  image. 

The  exact  length  of  the  frame  storage  buffer  a:  so  depends  on  the  assumed 
regions  of  support  for  the  impulse  responses  afn^, n2)  and  cfn^,  n2).  If  a(n1#  n2) 
and  cin^,  n2)  are  3x3  with  first  quadrant  support,  then  the  length  of  the  frame 
buffer  is  equal  to  the  number  of  pixels  ir*.  an  image.  However,  if  afn^,  n2)  and 
cin^,  n2)  are  3x3  but  centered  on  the  origin,  the  frame  buffer  length  must  be 
shortened  by  one  row  plus  one  pixel  to  ensure  that  the  forward  pixel  stream 
and  the  feedback  pixel  stream  are  properly  synchronized.  This  detail  further 
reduces  the  flexibility  of  this  structure  for  implementing  a  variety  of  iterative 
filters  on  a  variety  of  image  sizes. 

5.3  Row  Parallel  Architectures 

The  structure  for  implementing  an  iterative  filter  with  a  row  parallel 
image  flow  is  quite  simila.  in  style  tc  the  sequential  pixel  structure  discussed 
above.  However,  there  are  some  differences  in  implementation.  For  example, 
there  is  no  longer  a  need  for  the  shift  register  delays  in  Fig.  21  since  input 
pixels  are  available  in  parallel.  Also,  the  generation  of  the  output  pixels  in 
parallel  requires  the  replication  of  the  coefficient  multipliers  and  the  accumu¬ 
lator  in  Fig.  21,  Thus,  the  corresponding  structure  for  the  implementation  of 
an  FIR  filter  with  rev  parallel  image  flow  takes  the  form  shown  In  Fig.  22. 

The  brackets  represent  the  multiplication  by  filter  coefficients,  gating,  and 
accumulation  functions  shown  explicitly  in  Fig.  21.  The  brackets  also  serve 
to  indicate  the  inter- \x>w  communication  which  is  necessary  to  implement  a 
finite -extent  convolution. 
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The  structure  of  Fig.  22  may  be  substituted  for  the  boxes  labeled  Afw^, 
and  C(Wj,  in  Fig.  20  to  implement  an  iterative  filter.  Because  of  the  row 
parallel  image  flow,  the  number  of  data  paths  in  the  iterative  filter  is  equal  to 
the  number  of  rows.  This  number  may  be  prohibitively  large  for  many  applica¬ 
tions;  pixels  may  have  to  be  multiplexed  to  share  a  smaller  number  of  data 
paths.  This,  of  course,  leads  us  back  in  the  direction  of  the  sequential  pixel 
image  flow. 

The  row  parallel  structure  also  suffers  from  the  inflexibilities  described 
in  Sec.  5.2.  The  number  of  unit  delays,  the  amount  of  inter-row  communica¬ 
tions,  and  the  number  of  coefficient  multipliers  depend  on  the  size  of  the  con¬ 
volutions  being  implemented.  The  size  of  the  frame  buffer  in  Fig.  20  depends 
on  the  image  size  as  before. 

We  shall  just  mention  the  frame-by-frame  image  flow  in  passing,  since 
the  number  of  parallel  data  paths  required,  equal  to  the  number  of  pixels  in 
an  image,  is  impossibly  large  for  current  technology.  Each  output  pixel  in  a 
given  frame  is  computed  in  parallel  with  all  the  others.  For  each  convolution, 
this  requires  a  set  of  coefficient  multipliers  and  accumulators  for  each  pixel. 
Reducing  the  number  of  multipliers  and  accumulators  by  multiplexing  essen¬ 
tially  leads  us  back  toward  the  row  parallel  and  sequential  pixel  image  flows. 

5,4  Multiprocessor  Architectures  and  Algorithm  Partitioning 

The  architectures  discussed  in  Sec.  5,2  and  5.3  were  derived  from  the 
block  diagram  of  Fig.  20,  They  seem  to  possess  a  certain  inflexibility,  or 
hardwired  characteristic,  that  potentially  limits  their  usefulness,  although  a 
clever  designer  may  be  able  to  introduce  more  flexibility  into  the  structure  by 
using  random  access  buffers  and  microprocessor-controlled  multiplexing  of 
the  multiplication  and  accumulation  hardware. 

Alternatively,  we  can  try  a  software -oriented  approach  to  the  problem  of 
implementing  an  iterative  filter.  A  single  processor  with  an  appropriate 
amount  of  random  access  memory  can  be  straightforwardly  programmed  to 
cycle  through  the  necessary  computations  to  produce  the  output  image  frames 
one  pixel  at  a  time.  Now  suppose  we  have  a  number  of  processors  at  our  dis¬ 
posal.  The  problem  becomes  one  of  partitioning  the  image  processing 
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algorithm,  in  this  case  the  implementation  of  an  iterative  filter,  so  that  each 
processor  is  kept  busy  and  is  working  as  independently  as  possible. 

Realistically,  the  number  of  processors  will  be  relatively  small,  say  8. 

We  can  partition  the  image  into  strips  or  blocks  and  assign  to  each  processor 
the  responsibility  for  computing  the  output  pixels  in  a  particular  block.  If  the 
size  of  the  block  is  large  compared  to  the  extent  of  the  FIR  impulse  responses 
used  in  the  iterative  filter,  each  processor  can  run  independently,  as  if  it  were 
a  single  processor  working  on  a  smaller  image,  to  compute  the  necessary 
convolutions  on  each  iteration.  Care  must  be  taken  where  the  blocks  abut; 
some  inter-processor  communication  is  necessary,  either  to  transfer  partially 
computed  output  pixels  or  extra  input  pixels  to  neighboring  processors. 

If  one  of  the  processors  breaks  down,  then  one  block  in  each  output  image 
frame  will  be  garbled  or  missing.  (This  block  could  be  moved  around  from  one 
frame  to  the  next  by  reassigning  processors.)  Consequently,  this  partitioning  of 
the  problem  does  not  lead  to  a  very  robust  implementation.  As  an  alternative, 
we  could  subsample  the  input  image  frame  to  obtain  several  smaller,  lower- 
resolution  subimages.  A  processor  could  be  assigned  to  process  each  sub¬ 
image,  and  the  full  output  image  could  be  constructed  by  appropriately  inter¬ 
lacing  the  processed  subimages.  For  example,  Fig. 23  shows  how  an  8X8  image 
may  be  partitioned  among  p  =  8  processors.  Now  each  processor  has  only  one- 
eighth  of  the  input  pixels  to  handle  but,  for  a  finite-extent  convolution,  it  must 
generate  a  term  to  contribute  to  each  output  pixel  in  the  full  image.  Each  term 
is  roughly  one-eighth  as  complicated  as  the  full  sum  of  products  needed  to  com¬ 
pute  a  single  output  pixel,  so  the  amount  of  computation  performed  by  each 


Fig.  23.  Partitioning  the  pixels 
in  an  image  among  8  processors. 
The  numbers  indicate  the  pro¬ 
cessor  (p#  to  p7)  to  which  each 
pixel  is  assigned. 
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processor  is  essentially  the  same  as  if  it  were  computing  an  output  subimage 
independently  of  the  other  processors. 

This  notion  can  be  made  a  little  more  concrete  with  the  example  of  a 
3x3  convolution.  The  output  pixel  at  (n^,  n2)  is  given  by 

w(n1,n2)  =  a(-l,-l)  x(n^  +  l,n2  +  2)  +  a(-l,0)  x(n^  +  l,n2) 

+  a(-l,  1)  x(n^  +  l,n2  -  i)  +  a(0,-l)  xfn^,^  +  1) 

+  *(0,0)  x(n1,n2)  +  a(0, 1)  x(n1>n2  -  1) 

+  a(i,-i)  x(n4  -  i,n2  +  i)  +  a(i,0)  x(n4  -  i,n2) 

+  a(l,  i)  x(n1  -  l,n2  -  1)  .  (32) 

If  the  input  image  is  partitioned  as  in  Fig.  23,  then  processor  pjf  must 
compute  the  following  terms  and  send  them  to  the  other  processors: 

Source  Destination 


Processor 

Processor 

Term 

For  Output  Pixel 

Pff 

P0 

a(0,0)  x(n1#n2) 

w(n4, n2) 

pi 

- 

- 

p2 

a(i,-l)  x(nj  -  2,n2  +  2) 

+  a(-i,  1)  x(nltn2) 

w(n4  -  l.nj  +  i) 

P3 

a(l,  1)  x(n4,  n2) 

+  a(-l,-i)  x(nA  +  2,n2  +  2) 

w(n^  +  l,n2  +  1) 

p4 

a(l,  0)  x(n4,n2) 

w(n^  +  l.ng) 

p5 

a(-l,  0)  x(n1#n2) 

w(nj  -  1,  n2) 

p6 

a(0,~i)  x(n1#  n2) 

w(n1,n2-  1) 

P? 

a(0, 1)  x(n1#n2) 

w(n4,n2  +  1) 
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Similarly,  processors  pi  to  p7  must  deliver  partial  results  to  all  8  proces¬ 
sors.  Then  the  partial  results  can  be  accumulated  to  generate  the  output  pixels 
assigned  to  each  processor. 

Let  us  look  briefly  at  the  amount  of  memory  required.  We  shall  assume 

2 

that  image  frames  consist  of  N  pixels  and  there  are  P  processors.  Each 

2  2 

processor  needs  an  N  /P-pixel  input  buffer  and  sin  N  /P-pixel  output  buffer. 

In  addition,  each  processor  must  have  N  (P  -  1)/P  storage  locations  for  saving 

the  partial  results  which  will  be  sent  to  the  other  (P  -  1 )  processors.  Conse- 

2 

quently,  the  total  amount  of  storage  required  for  all  processors  is  N  (P  +  1) 
pixels. 

When  the  processors  are  done  computing  the  partial  outputs,  they  need  to 
send  and  receive  data  to  allow  them  to  finish  computing  the  output  pixels.  We 
need  to  postulate  an  interprocessor  sv/itch  to  accomplish  this  communication. 
The  p0  output  is  coupled  to  the  pi  input  while  the  pi  output  is  coupled  to  the 
p2  input,  and  so  on.  Then  p0  communicates  with  p2,  pi  with  p3,  and  so  on. 
After  P  -  1  settings  of  the  switch,  the  final  output  pixels  can  be  computed.  The 

7 

total  number  of  pixels  transmitted  over  the  interprocessor  switch  is  N  (P  -  1), 
but  because  of  the  parallelism  only  hi  (P  -  i)/P  transfer  cycles  are  used. 

The  amount  of  memory  used  by  each  processor  to  buffer  the  intermediate 
results  can  be  reduced  if  the  intermediate  results  can  be  transmitted  shortly 
after  they  have  been  computed.  Obviously,  the  number  of  pixels  transferred 
and  the  number  of  transfer  cycles  needed  do  not  change;  more  frequent  trans¬ 
fers  of  shorter  data  blocks  occur. 

5.5  Conclusions 

The  ideas  outlined  herein  should  be  regarded  as  preliminary.  Neverthe¬ 
less,  some  conclusions  can  be  drawn.  First,  the  development  of  a  structure 
for  implementing  Iterative  filters  based  on  the  block  diagram  In  Fig.  20  seems 
to  lead  to  rather  inflexible  architectures.  Second,  for  particular  applications, 
the  type  of  image  flow  may  depend  on  the  type  of  sensor  (flying  spot  vs  linear 
array  vs  lens).  Multiplexors  and  demultiplexors  may  be  used  to  translate 
from  one  image  flow  to  another,  of  course,  but  processing  architectures  may 
have  to  be  sensitive  to  the  differences  in  image  flows  in  a  real-time  application. 
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At  the  moment,  there  seem  to  be  two  interesting  methods  for  partitioning 
the  filtering  operation  which  deserve  closer  scrutiny.  The  first  involves 
assigning  each  processor  to  compute  the  output  pixels  in  a  particular  block  of 
the  output  image  frame;  the  second  assigns  a  subsampled  image  to  each  pro¬ 
cessor  so  that  the  effects  of  a  malfunctioning  processor  will  be  more  diffuse 
and  perhaps  less  deleterious. 
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